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Abstract 

Planets are built from planetesimals: solids larger than a kilometer which grow by colliding pair- 
wise. Planetesimals themselves are unlikely to form by two-body collisions; sub-km objects have 
gravitational fields individually too weak, and electrostatic attraction is too feeble for growth 
beyond a few cm. We review the possibility that planetesimals form when self-gravity brings to- 
gether vast ensembles of small particles. Even when self-gravity is weak, aerodynamic processes 
can accumulate solids relative to gas, paving the way for gravitational collapse. Particles pile 
up as they drift radially inward. Gas turbulence stirs particles, but can also seed collapse by 
clumping them. While the feedback of solids on gas triggers vertical shear instabilities that ob- 
struct self-gravity, this same feedback triggers streaming instabilities that strongly concentrate 
particles. Numerical simulations find that solids ~10-100 cm in size gravitationally collapse in 
turbulent disks. We outline areas for progress, including the possibility that still smaller objects 
self-gravitate. 
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1 INTRODUCTION 

The nebular hypothesis — that planets coalesce from disks of gas and dust orbiting 
young stars — is confirmed today in broad ou tline. Protoplanetar y disks of ages 
1-10 Myr are now studied routinely (Figure [T] I Watson et al. 2007 ). They contain 



enoug h mass to spawn planets like those detected around hundreds of Gyr-old 
stars ( Jones et al. 20081 ). Debris disks, observable via dust generated from colli- 



sions bet ween larger p arent bodies, bridge our understanding in the 10-100 Myr 
interval dWvatt jg . These parent bodies may represent first-generation plan- 



etesimals, the building blocks of planets. Our Solar Syste m may preserve a record 
of planet formation, in t he size distributio ns of asteroids ( Morbidelli et al. 20091 ) 
and Kuiper belt objects ( Pan &; Sari 20051 ). 



The remarkable journey that disk solids make in growing from microscopic 
dust to Earth-mass (-M©) planets divides into three legs. At the smallest sizes 
<J cm, chemical bonds and van der Waals forces enable grains to stick to one 
another. At the largest sizes 3> km, gravity promotes growth. Pairs of objects 
coagulate upon colliding because their gravity is strong enough to retain collision 
fragments. 

At intermediate sizes lies the domain of planetesimal formation. It presents 
the most challenging of terrains; here gr ain surfaces are insufficiently sticky and 
the gravity of a single object too feeble. ISafronov kindled the hope that 



where the individual fails, the collective might succeed: if solid particles have a 
large enough density en masse, their gravity can draw them together. The many 
calculations seeking to realize this vision are the subject of our review. 

An incompressible fluid body in hydrostatic equilibrium, orbiting a star of 
mass M* at distance r, has sufficient self-gravity to resist disruption by the 
star's tidal gravitational field if its density exceeds the Roche value /0R oc h c ~ 
3.5M*/r 3 (jChandrasekhar 19871 ) . Though we will see that the effects of self- 



gravity manifest at arbitrarily low densities ( §5.2)) . the Roche density is a fine 
benchmark against which to measure progress towards forming planetesimals, 
and we will use it as such. It is a formidably large density, 2-3 orders of m agni- 
tude greater t han typical densities at di sk midplanes (32J). Safronov (19691 ). and 



independently Goldreich Ward (19731 . hereafter GW), proposed that densities 



approaching Roche (actually a factor of 20 smaller, see §5.ip could be achieved 
by hav ing dust settle vertical ly and accumulate in a thin "sublayer" at the mid- 
plane. IWeidenschilling pointed out that large densities were not so easily 



achieved — that turbulence generated by vertical shear across the sublayer would 
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halt settling and prevent gravitational instability (GI). GW appreciated that the 
sublayer would be turbulent, but overlooked the decisive impact of turbulence on 
the layer's thickness and density. 

This stalemate has been broken — indeed the entire playing field redrawn — in 
the last decade by abandoning three assumptions. First, the height-integrated 
surface density ratio of dust to g as need not be solar. Enriching the disk in metals , 
say by photoevapora tion of gas (jThroop Bally 2005, iGorti Hollenbach 2009, 
Ercolano et al. 20091 ) or by radial drifts of particles ( Youdin Shu 2002, here- 
afte r YS), can st op vertical shear turbulence from forestalling GI ( Sekiva 19981 . 
YS, Chiang 20081 ). Second, turbulence, of whatever origin, is not always the en- 
emy of GI. Magneto-rotational tubulence is observed in nu merical simulations to 
generate long-lived str uctures in gas that can trap particles (jFromang fc Nelson 20051 . 
Johansen et al. 2006bl ). Third, gas drags particles, but by Newton's Third Law, 
particles also drag gas. Properly accounting for the backreaction of particles on 
gas leads to powerful drag i nstabilities that can concentrate particles and seed G I 
(Goodman fc Pindor 2000 . lYoudin fc Goodman 200fil . I.Tohansen fc Youdin 20071 ) . 



This review is organized as follows. Section 2 establishes order-of-magnitude 
properties of protoplanetary disks using an up-to-date model for the minimum- 
mass solar nebula. Attention is paid to primitive chondritic meteorites. Section 
3 gives a primer on the aerodynamics of grains, treating gas drag in the "test 
particle" limit where backreaction is neglected. We discuss how grain dynamics 
are affected by gas turbulence. Section 4 describes the extent to which grains 
can grow by sticking. Section 5 rehearses and comments on Toomre's criterion 
for GI — and how the criterion is removed when self-gravity is combined with gas 
drag. Section 6 treats sublayer shearing instabilities. Section 7 introduces secular 
drag instabilities that concentrate particles w ithou t recourse to self-gravity. We 
offer new insights into Goodman & Pindor 's ( 2000l ) toy model for drag instabil- 
ities. Section 8 outlines the nonlinear outcome of GI, describing breakthrough 
numerical simulations. Finally, a summary is supplied in Section 9, where we also 
list a few forefront problems. Our study is peppered throughout with calculations 
that readers are encouraged to reproduce — and improve. 

Our review complements others in the field of planet formation. Blum fe Wurm (20081 ) 
do more j ustice than we have to the growth of particle aggregates by grain-grain 
sticking. iGoldreich et al. f2004h provide a pedagogical and cutting-edge review 
of protoplanet accretion by gravitationally focussed, pairwise collisions. For a 
review of planetesimal formation that overlap s ours and co vers some topics in 
greater detail, see the Les Houches lectures by lYoudin (20081) . 



2 PROTOPLANETARY DISKS 

A starting point for calculations is the minimum- mass solar nebula (MMSN), de- 
rived by adding enough H and He to solar system planets to restore them to solar 
composition, and spr eading the augmented m asses into abutting annuli centered 
on their orbits (e.g., Weidenschilling 1977bl ). This exercise is uncertain, since 



estimates of the total "metal" (non-H and non-He) content of Jup iter range from 



upit 

10 to 42 M e ; for Saturn the range is 15 to 30 M ffi (jGuillot 20051 ) . Even if these 



uncertainties were reduced, it is unlikely our solar system formed strictly from 
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the MMSN. The conversion of disk metals into planets cannot be 100% efficient, 
and current pla netary orbits may differ significantly from their original ones (e.g., 
Malhotra 1993 ). Moreover a major theme of this review is that planetesimal for- 



mation may require the metal fraction to evolve to values greater than that of 
the bulk Sun, at least in certain regions of the disk. Thus the utility of the 
MMSN lies in providing a baseline for discussing real-world complications, and 
in establishing orders of magnitude. 

For numerical estimates in this review, we adopt a disk surface density: 

/ \ -3/2 

S g = 2200F(^— J gem- 2 (1) 
/ r \-3/2 

S p = 33FZ rcl ^— J gem" 2 (2) 

where subscripts g and p denote gas and particles (condensed metals), respec- 
tively. Our MMSN [F = 1, Z ve \ = 1), uses the updated condensate mass fraction 
for Solar abundances of 0.015 ( Lodders 20 03). When methane ice sublimates 



above temperature T ~ 41 K, Z VP \ = 0.78; abo ve ~182K, water and all other 
ices are lost so that Z re \ = 0.33 (|Lodders 20031 1. The coefficient for eq. ([2]) is 



chosen to give 1M® of solids in an annulus centered on the Earth's orbit when 
F = 1 and Z re \ = 0.33. We invoke values of Z re \ > 1 to account for various metal 
enric hment proces ses, depicted in Figure [2] an d discussed throughout this review 



fe.g.. lSekiva 19981 . YS. lThroop fc Bally 2 005) 



Integrated to r = 100 AU, eq. ([I]) yields a mass of 0.03FM©. Astronomical 
observations of disks orbiting Myr-old, Sun-like stars suggest they contain 0.001- 
0.1 M Q = 1-100 Mj of gas and dust (| Andrews fc Williams 20051 ). Disk masses 
are derived from mm-wave radiation from dust, located ~100 AU from host 
stars. From the dust emission is calculated a dust mass, and from the dust mass 
a gas mass is extrapolated assuming a solar dust-to-gas ratio. The gas mass thus 
imputed is uncertain because the disk metallicity is not known, and because the 
modeled dust mass depends on an unknown grain size distribution. Disk masses 
may be systematically underestimated because grains may have grown to sizes 
S> mm and would the refore be practically invisible at sub-cm wavelengths (e.g., 
Hartmann et al. 20061 ). 



Based on their near-ultraviolet excess emission ( Calvet Sz Gullbring 19981 ). youn. 



Sun-like stars accrete gas from disks at a typical rate ~ 10 M Q yr (e.g., 
Hartmann et al. 20061 ). Accretion implies that disk gas cannot everywhere be 



static ( "passive" ) . Mass is transported inward by the ou tward transport of angu- 
lar momentum, either by turbulence or by ordered flows ( Frank et al. 20 02). The 



transport mechanism — or, as is commonly stated, the origin of disk viscosity — 
remains obscure, with prop osals ranging from turbule nce driven by th e magneto- 
rotational inst ability (MRI; Balbus 20091 ') to vortices ( Lithwick 20091 ) to gravita- 



tional torques ( Vorobyov &: Basu 20081 ). When transport is local, our ignorance 



is encapsulated in the parameter a (e.g. JFrank et al. 20021 ). the ratio of the local 
shear stress to the total pressure. Crude, disk-averaged val ues of a ~ 10 -2 are 
inferred from observations of M* vs. age (jCafvet et al. 200(1 ). 



Accretion may be restricted to certain radii and times, as in models that rely 
upon the MRI — not all of the disk may be suffici ently ionized to couple to mag - 



netic fields (e.g., iBai &: Goodman 20091 . see also IChiang Sz Murrav-Clav 200 
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Magnetically inert regions are called "dead zones" ( Gammie 19961 ). Whether 



planetesimals form in active or passive regions of the disk is an outstanding ques- 
tion. This review will explore both cases. 

The fraction of stars with near-infrared excess emission attributable to optically 
thick, gase ous disks decreases from near unity at stell ar ages <^ 1 Myr, to 5% at 
£ 10 Myr (jHernandez et al. 20081 . Iffillenbrand 20051 ). Whatever combination of 



accretion, photoevaporation, and planet formation is responsible for this observed 
decline, giant planets — including ice giants, whic h contain more hydrogen th an 
can be explained by accretion of hydrated solids ( Lissauer Stevenson 20071 ) — 



must form within several Myr before the gas dissipates. The same deadline 
characterizes planetesimal formation, since the ice giants, which by mass are 
80-90% rock and ice, coagulated from planetesimals. 

Current models set the formation of solar system giant planets within a disk 
with more solids, and sometimes more gas, than the MMSN. lGoldreich et al. f2004T ) 



identify Neptune and Uranus with "isolatio n-mass" oligarchs at r ~ 20-30 AU, 
requiring FZ TC \ ~ 60 iLissauer et al. f2009h prefer FZ TC \ ~ 3.5 at r ~ 5AU, so 
that Jupiter's core can accrete its gaseous envelope within 3 Myr. 

At the midplane of a passive disk — one heated solely by stellar radiation — the 
gas temperature, scale height, and density are approximately 

/ r \-3/7 

r = 120 Uu) K (3) 

/ r \2/7 

K = 0-022r (-) (4) 

/ r \ -39/14 

p g = 2.7 x 1(T 9 F ^— J gem" 3 . (5) 



These are adapted from lChiang &: Goldreich (19971 ). adjusted for a disk obeying 



CD - ©, orbiting a young star of mass M* = 1M , radius i?* = 1.7 R & , and 
temperature T* = 4350 K. Thus in our MMSN, water ice starts to condense 
outside ~0.4 AU, and methane f reezes outside ~12 AU. More detailed models 
of passive disks are reviewed by iDullemond et al. (20071 ). Turbulent accretion 



can give higher midplane temperatures. The density and thermal structure of 
ac tive disks depends on the assumed viscosity profile, a(r). One model is that 



of iD'Alessio et al. (20011 ). which assumes a{r) = constant and is tailored to fit 
broadband spectra and reflected light images. For simplicity, estimates in this 
review employ (pQ)-([5]) for both active and passive disks. 

What do astronomical observations tell us about grain sizes? In disk surface 
layers directly illuminated by optical light from host stars, grain sizes are 1-1 /tm, 
as deduced from mid-infrared silicate emission b ands dNatta et al. 20071 ). and 



from scattered light images at similar wavelength ( McCabe et al. 2 003). Surface 
grains have settled vertically, residing at heights z above the midplane of l-3h g 
(jChiang et al. 200lh : :'or a disk in which dust and gas are well mixed, z/h g ~ 4-5. 



At disk midplanes, mm to cm-sized grains are routinely invoked to match mm to 



A protoplanet of mass M v becomes "isolated" when it has accreted an annulus of disk 
material having a wid th about 5x the radius of its Hill sphere, -Rh ~ [Mp/fSM,)] 1 '^ ( e -g-, 
IGreenberg et al. 199ll ). The average density of the planet spread through its Hill sphere is of 
order pR OC h G ~ 3.5M,/r 3 . 
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cm-wave spectra and images (e.g., iD'Alessio et al. 200ll . iTesti et al. 20031 ) 



The most primitive meteorites, of nearly solar photospheric composition except 
in volatile elements, also offer data on particle sizes. Confounding meteoriticists 
and astronomers alike is why up to ~90% of their volume is filled with chon- 
drule s: once molte n, 0.1-mm to cm-sized spheres that solidified 4.57 Gyr ago 
(e.g., Hewins 19961 ). Chondrule petrology is consistent with their having been 



heated just above liquidus for less than minutes, and having cooled for hours to 
days. The heating mechanism is not known; nebular s hocks are suspected, but 
the origin of such shocks is debated (|Desch et al. 20051 ). Chondrules might well 



have been the building blocks of the first-generation planetesimals, brought to- 
gether by self-gravity. That ~10% of chondrules are binaries which collided and 
fused while still partially molten implies th at when chondrules w ere suspended 



in space, they had large collective densities (jGooding &: Keil 19811 ). possibly ex- 
ceeding the Roche value. Extremely d usty environments are also indicated by 
the retention of volatiles in chondrules ( Alexander et al. 20081 ). 



3 AERODYNAMICS OF INDIVIDUAL PARTICLES 

Gas drags particles. The degree of coupling is measured by the dimensionless 
stopping time r s = SIk^s, where f?K is the Keplerian angular velocity and 

t = mv j Fn ~ J Pss/(p g c g ) if s £ 9A/4 (Epstein) 

s - mv ^'^ ~ \ 4p sS 2 /(9p gCg X) if s % 9A/4, Re £ 1 (Stokes) W 

measures how long it takes a particle of mass m, radius s, and internal den- 
sity p s to have its speed v re \ relative to gas be reduced by order unity. The 
gas sound speed is c g . Particles are well entrained in gas when r s <C 1. We 
have given the two cases for th e drag force -Frj most relevant for planetesi- 
mal formation: Epstein's ( 19241 ) law of free molecular drag, arising from the 



difference ~/3 g [(c g + v ve \) 2 — (c g — w re i) 2 ] i n momentum fluxes received by the 
windward and leeward faces of the particle, and Stokes drag for low Reynolds 
number Re = st> re i/(Ac g ). Given ([5]), the mean free path for collisions between 
gas molecules is A 0.5 F~ 1 (r/AU) 39 / 14 cmH O ther cases for Fd are given by 
Adachi et al. (19761 ) and IWeidenschilling (1977al ). For F = 2 and p s = 1 g/cm 3 , 



marginally coupled bodies (r s = 1) have s = s± increasing from 35 to 120 cm in 
the Stokes regime as r runs from 1 to 7 AU; thereafter s± declines with r in the 
Epstein regime, with s\ ~ 13 cm at 30 AU. 

Gas and dust move at different velocities because pressure gradients barely 
accelerate particles with p s 3> p g . For the moment, let us neglect turbulence, and 
assume the collective particle density p p -C p g so we can ignore the backreaction 
of dust on gas. Then the azimuthal gas velocity v g s obeys 



vL GM* 1 OP 



r r 2 p g dr 



(7) 



2 We use a constant molecular cross section, cr(H2) ~ 2 x 10 15 cm 2 , that reproduces the 
dynamical viscosity fi at 200 K. But constant ct(H2) assumes /j, oc VT, which fails for T < 70 
K, the Sutherland constant for H2. At colder temperatures ji oc T 3 / 2 so that <r(H2) oc T _1 
l|Chapman fc Cowling 19701 ) . For simplicity we neglect the lower values of viscosity and A that 
this effect would produce in cold regions. 
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Pressure P ~ PgC 2 . makes the gas rotate more slowly than the local Keplerian 
velocity vk = ^k?~ = \J GM*/r by 

dP/dlnr ( r X 1 / 14 
ijvk = ^g<A ~ ^ ~ 25 — m s (8) 

using the thin disk approximation c g /i>K = ^gA ^ 1- While t/vk ~ c 2 /vk is 
independent of total disk mass and is nearly independent of r, it increases in 
hotter disks, perhaps to ^ 50 ms -1 . An individual particle seeks to orbit at the 
full vk- Thus it normally experiences a headwind, except possibly where there 
are variations to the power-law descent of P(r), as considered at the end of §3.1.21 



3.1 Particle Drifts and Consequences 

3.1.1 Laminar Drift Speeds In a passive disk, a particle settles vertically 
to the midplane and, on longer timescales, drifts radially i nward as the h eadwind 
saps its angular mom entum. The following is drawn from lYoudin (20081 . see also 



Nakagawa et al. 19 86). If we neglect backreaction, the cylindrical components of 



the gas velocity are v gr = 0, v gz = 0, and v g< f, = (1 — r^Sl^r- In z, the particle's 
equation of motion reads 

z = -z/t s - Q K z (9) 

where the last term accounts for stellar gravity, for z <C r. For r s <C 1, the 
particle settles at terminal velocity — Q K zt s . For r s 3> 1, the particle behaves 
as a lightly damped harmonic oscillator whose amplitude decays as e~ t ^ 2ts \ A 
characteristic time for settling, valid for all r s , is 



t '^{-^)- (10> 

In r and (ft, 

r — rift = —v K /r — r/t s (11) 

r4> + 2r<j) = -(r<i)-v g4> )/t s . (12) 

Writing (ft = Ok + Sv^/r where \6v^,\ <C Ok?", we approximate (ft Ok ~ 
— 3Sli<T/(2r) and drop \r\ <C |^As|- Both approximations can be checked a pos- 
teriori. To first order in 



6v^ = {4>- n K )r « . (14) 

The time for radial drift, t r = \r/r\, is longer than t z . Both times are 

minimized, and velocities relative to gas are maximized at ~??t>K; for marginally 
coupled bodies. This is the well-known problem that boulders having sizes si ~ m 
drift towards the star in min(t r ) ~ (t^Ok) -1 ~ 200(r/ AU) 13//14 yr, too quickly to 
form planets. 
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3.1.2 Pileups AND Pressure Traps As par ticles of a given size dr ift in- 
ward, they tend to "pile up," increasing S p /S g (YS. lYoudin Chiang 20041 ); see 
Figure [2j Their inward mass flux decreases with decreasing r if Epstein drag ap- 
plies and the outer regions are not already depleted in solids. Idealized pile-ups 
march inward on the drift timescale with increasing amplitude. The actual abil- 
ity of particles to pile up coherently will be affected by ongoing evolution of the 
particle size distribution and by any disk turbulence. Youdin Chiang (20041 ) 
explore one prescription for turbulent viscosity predicate d on the KHI that am - 
plifies pile-ups, but its bas is is unproven; see §5.2.3 of iGaraud &: Lin (20041 ) . 
Stepinski &: Valageas (19961 ) show, by combining the drift and turbulent diffu- 
sion of particles (see £ )3.2p with the viscous evolution of £ g , that the solids-to-gas 
ratio in a disk evolves, and often increases, as it accretes onto the star. 



A special pile- up of dust could occur just inside an ice- line (e.g. Ida &: Lin 20081 ). 
When "dirty snowballs" drift towards and evaporate i nside the ice-line, they un- 
leash small dust grains that may accumulate there fsee lCiesla fc Cuzzi 20061 . who 
included this effect in a coagulation model). Furthermore, the ice-line may act 
as a "cold trap": vapor t hat diffuses radially outwa rd and crosses the line will 
condense and accumulate ( Stevenson &; Lunine 19881 ) . Though water is the dom- 
inant volatile, the methane condensation front could also be significant (32]). 

We have shown that inward particle drift is a direct consequence of gas pres- 
sure decreasing radially outward. By the sa me physics, sol ids of all sizes drift into 
and collect within local pressure maxima (jWhipple 19721 ). with r s ~ 1 particles 
accumulating fastest. If a pressure bump has a narrow radial width I < r, even 
small amplitudes 5P ~ (£/r)P will produce local maxima. P roposed sources of 
press ure bumps inclu de gaseous spiral arms (jRice et al. 20041 ). anticyclonic vor- 
tices ( Chavanis 2000l ). the aforementioned ice - lines, and fluctuations in magneto- 
rotational turbulence (IJohansen et al. 2009a. sec £|3.2.3p . including spatial varia- 
tions in the backgro und flux ( Kato et al. 20091) or m agnetic resistivity, as at the 
edges of dead zones (|Gammie 19961 ) lYoudin (20081 ) discusses pressure trapping 
in more detail. 



3.2 Turbulent Stirring 



The response of solids to turbulent gas is multi- faceted (see lToschi k, Bodenschatz 20091 . 
especially their section 4.1). Over many eddy times, turbulence diffuses solids, 
vertically and radially. On shorter timescales, transient flow structures can have 
just the opposite effect: they can collect particles. 

Many studies assume that turbulence in disks is Kolmogorov in character, and 
in our example calculations below we also adopt this view. We associate 5v , 
t Q , and i ~ 5v t with the speed, turnover time, and length of the largest and 
fastest "outer scale" eddies. Smaller eddies of length scale i have slower speeds 
5v£ ~ <5f (£/^ ) 1//3 and shorter turnover times te = £/5v£ tx £ 2 / 3 . The cascade of 
energy from larger to smaller scales terminates at the "inner scale," characterized 
by l\ ~ v^/ A t l J 4: /5v^ 2 , t{ ~ yjvt /5v , and 5v\ ~ l\jt\ set by the molecular 
viscosity v. For many estimates we assume that outer scale eddies have lifetimes 
limited by orbital shear, so that t Q ~ Sl^ 1 . By ignoring the difference between the 
turbulent diffusivity for mass, D g ~ 5v £ ~ 5v^t , with the turbulent diffusivity 
for angular momentum, ac g h g ~ aCgfi^ 1 , we extract the convenient prescriptions 
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5v Q ~ -y/a-Cg and £ Q ~ y/ah % ( Cuzzi et al. 200ll . see also lYoudin Chiang 200 . 



However the assumption of homogeneous, isotropic turbulence characterized 
by cascades of energy down to smaller scales may be i ncorrect. Simulations 
of MR I turbulence, both glob al ( Fromang &; Nelson 20051 ) and in large shearing 



boxes ( Johansen et al. 2009a! ) . show that inverse cascades produce pressure per- 



turbations that survive many tens of orbits and have in-plane sizes larger than 
h g . These coherent structures may accumulate particles by pressure trapping 
(Figured 33X21 . 

3.2.1 Particle-Particle Velocities In active regions of the disk, rela- 
tive gas-particle and particle-partic le velocities may be dominated by turbulence 
dKusaka et al. 1970l.lCameron 19731 1. rather than the systematic drifts considered 



m 



33XT1 IVoelk et al. (1 980) derive integral expressions for these relative ve- 



locites, evaluating them numerically for Kolmogorov turbulence. Consider iden- 
tical grains with t{ < t B < t Q , i.e. stopping times between the inner and outer 
scales. Collision speeds are acquired from eddies with a turnover time tg ~ t s . If 
t£ -C i S i coupling to the eddy is too weak to excite particle motion. If t£ t s , 
grains are so tightly coupled to the incompressible gas flow that they do not 
collide. Thus colliding grains have relative velocity 5v P)P ~ 5v£ ~ 5v {t s /t ) l l 2 — 
which is not to be confused with the random speed of an individual particle, 
5v p ~ 5v Q . 

For t s > t Q , particles are loosely coupled to all eddies. During t s there are 
t s /t Q kicks from the strongest, outer scale eddies, so the random wa lk in particle 
velocity saturates at 5v p ~ 5v Q / \Jt s /t Q ( Youdin <fe Lithwick 20071 ). Since the 



random walks of different particles are uncorrelated, the random velocity is now 
also the typical relative speed. 

Combining our expressions for collision speeds gives 



H,p ~ <Wst/(i + st 2 ) (15) 



where the Stokes numb er St = t s /t Q . This agrees roug hly with lVoelk et al. (19801 ) . 



Weidenschilling (19841 ) provides more precise fitting formulae — further refined by 
Ormel Cuzzi (2007) — that include important corrections for t s < t\ and un- 



equal particle sizes. lYoudin &: Lithwick (20071 ) confirm that orbital dynamics 



introduces only modest corrections when t Q <^ K . 

Thus collision speeds in disks peak at 5v Q for r s ~ 1 particles, since St ~ r s when 
t Q ~ fi^ 1 . Magneto-rotational turbulence with 5v a ~ ^facg ~ 70(a/10 -2 ) 1 / 2 
(r/AU) _3//14 ms" 1 can induce collision speeds faster than drift speeds (eqs. [TBI 
and[U 



3.2.2 Diffusion and Particl e Scale Height On t imescales > t Q , tur- 



bulence diffuses particles spatially. lYoudin &: Lithwick (20071 ) compute, formally 
and by order-of-magnitude methods, the particle diffusivity D p as a function of 
the gas mass diffusivity D g . When r s is small, orbital dynamics are negligible. 
If further the Stokes number St = t s /t Q < 1, then particles are well coupled to 
gas and D p ~ D g . If St > 1, it takes t s for the loosely coupled particle to have 
its random velocity 5v p (derived in the paragraph just above eq. [T5l) changed by 
order unity. Then D p ~ 5v p t s ~ 5v^t ~ D g . Thus when t s ^ 1, D p ~ D g 
regardless of St. In general, D p ~ D g /(1 + r 2 /4) because orbital epicycles limit 



10 



CHIANG & YOUDIN— SUBMITTED DRAFT 



diffusion (jYoudin k Lithwick 200% . 



The time for particles to diffuse radially r 2 /D p is shorter than the drift time t r 
when a > t s . Thus for a ~ 10~ 2 , sub-cm-sized pa rticles are well coupled to the 
disk accretion flow; but see Takeuchi &: Lin (20021 ) for possible height-dependent 
complications. 

By equating the vertical diffusion time hp/D p with the gravitational settling 
time t z , we find a turbulent dust scale height 



"WiSr^ (16) 



for all t s . ICarballido et al. f2006) derive and numerically confirm ([TBI) for t s 3> 1 



The t s C1 limit is well-known ( Cuzzi et al. 19931 . Dubrulle et al. 19951 ). but see 



the end of §3.2.31 for a possible correction. 

3.2.3 Turbulent Concentration Between Eddies On timescales < 
t Q , particles tend to be centrifugally flung o ut of high vo r ticity eddies, and fun- 



neled into their interstices o f lower vorticity (jMaxey 19871 . lEaton &; Fessler 1994 . 



Toschi fc Bodenschatz 2009h . Since centrifugal support in eddy vortices makes 



them low pressure centers, turbulent concentr ation is consi stent with the ten- 
dency of particles to seek high pressure ( §3.1.21 Youdin 20081 ). Particles of given 



t s are concentrated preferentially by eddies that turn over on the same timescale; 
in other words, particles are concentrated most strongly by eddies to which they 
are marginally coupled. In a Kolmogorov cascade, smaller eddies are more ef- 
fective than larger eddies at concentrating their respective, marginally coupled 
particles. This is because vorticity ~ 6vi/£ oc l~ 2 ^ , down to the dissipation scale 
where h ~ 5(r/AU) 135 / 56 m and U ~ 30(r/AU) 9 / 4 s for F = 2, a ~ 0.01, and 
t D ~ ^k 1. Turbulent concentration is special since it can collect r s <C 1 particles. 
Those particles that are optimally concentrated are marginally coupled to the 
smallest eddies: they have t s = t\, which corresponds to s = s\ ~ 0.1 mm at 1 AU 
and Si ~ Wpm at 30 AU, with some uncertainty in the viscosity (see footnote [2]) . 

Cuzzi et al. (20oil ) apply turbulent concentration to chondrules: objects hav- 



ing sizes 0.1-1 mm are optimally concentrated at r = 2.5 AU in their model disk, 
consistent with our estimate of Si above. The shape of the chondrule size distri- 
bution is well reproduced by turbulent concentration. This raises two interesting 
possibilities: (1) the chondrule size distribution was initially much wider, and 
has since been narrowed to the one observed today by turbulent concentration; 
or less likely (2) chondrules happen to form with sizes that closely match optimal 
concentration scales at the location of the asteroid belt. 

The upper limit on concentration is set by mass loading, the backreaction of 
particles on th e turbulence. In direct numerical simulations of si-sized particles in 
turbulent gas, Hogan Cuzzi (20071 ) find a maximum concentration factor <3? = 



Pp/Pg ~ 100 at the grid, i.e. dissipation scale. This maximum value obtains for 
moderate Reynolds numbers Re Q = i 5v /u and a box-averaged ($) = 1. Extrap- 
olating these results using "cascade multipliers" ( Sreenivasan &; Stolovitzkv 19951 ). 



they show that max<I> ~ 100 still applies at higher Re Q . Maximal clumping at 
the dissipation scale yields a mass Qpgif equal to that contained in a compact 
10-cm solid at r = 2.5 AU. However these clumps are unlikely to compact- 
ify because of their extremely short characteristic lifetimes, on the order of i; 
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(|Eaton k Fessler 19941 1. For our disk parameters, i; ~ 4 minutes. 



Cuzzi et al. (20081 . hereafter CHS08) introduce the possibility that Si-sized solids 



concentrate on much larger length scales I* ~ 10 4 km, or 10 5 ^; at r = 2.5 AU. The 
clump mass ~$yOg^ would correspond to a compact ~20-km planetesimal, for 
~ 100 concentration. Large scale fluctuations might arise as an "intermittent" 
phenomen on: deviatio ns from self-similarity that show up in high order structure 
functions ( Frisch 19 96). Evidence exists for particle clumping at larger scales, i.e. 



in the inertial range of isotropic turb ulence, but preferentially for particles larger 
than Sj, with t s > t; ( Bee et al. 20071 ). This is to be expected since eddies in the 



inertial range have turnover times longer than t{. Whether clumping of Sj-sized 
or larger particles happens at astrophysically interesting amplitudes and rates is 
an active area of research. 

CHS08 find that an ^-scale clump is massive enough to survive ram pressure 
stripping as it plows through gas. They propose that self-gravity draws solids to 
the clump center at the terminal velocity, a process that lasts ~100 orbits. The 
turbulence that is invoked to create the clump could still destroy it, since the 
lifetime of an 4-eddy at 2.5 AU is i* ~ i (4/4>) 2/3 ~ 0.007 orbit. Streaming 
instabilities ( §7.2h may aid in clump survival, as they show that mass loading 
can promote particle clumping in Keplerian disks. The enticing possibility that 
chondrules comprise first-generation planetesimals warrants further investigation 
in this area. 

Eq. (|16p for the turbulent dust scale height might need rev ision because turbu- 
lent concentration accelerates vertical settling ( Maxev 19871 ). The correction has 



not been investigated, but would apply for St <C 1. Turbulent concentration is 
also neglected in studies of grain growth by sticking, a subject to which we now 
turn. 



4 PARTICLE GROWTH BY STICKING 

Particles stick upon colliding if they move slowly enough, dissipate enough en- 
ergy during impact, and are small enough — since surface area increases relative 
to mass for smaller bodies. We can estimate the relevant orders of magnitude 
by comparing the initial kinetic energy of two identical elastic spheres collid- 
ing at speed v co \, with the surface binding energy at the moment of maximum 
deformation. From Hertz's law of contact, the radius of the compressed cap is 
b ~ s(p s v^ ol / 'E) 1 / 5 , where E is Young's modulus. The binding energy is ~7& 2 , 
where 7 is the surface tension from unsaturated bonds, and exceeds the collisional 
energy for 



7 \ 5 / 6 / 7 x 10 10 erg cm" 3 \ 1/3 

^col < ^stick ~ 2.8 ——L ^ ^ 17 

370 erg cm z 1 



1 g cm 

Ps 




for ice (jChokshi et al. 199irYo"udin 20041 ). For jum-sized spheres made of silicate, 



(fT?]) yields v s tick ~ 7 cms , about an order of magnitud e smaller than sticking 
velocities measured experimentally ( Blum &; Wurm 20081 ). 
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In laminar disks, growth by sticking is expected to stall at ~cm sizes. From 
equation (fl~3j) . radial drift velocities exceed ~lms _1 — the sticking velocity for 
/im-sized monomers — once r s ~ 0.02, corresponding to particle sizes s ~ 5 cm, 
4 cm, and 2 mm at 1, 5, and 30 AU, respectively, for F = 2. Thus super- 
cm particles would fail to accrete smaller grains that move with the gas. As 
judged by equation (fT7|) . comparably-sized, super-cm bodies would stick only at 
unrealistically low relative speeds, with extremely low surface binding energies 
per unit mass. Indeed bodies >, m m in size have not been observed to stick 
experimentally ( Blum &: Wurm 20081 ) . 



Parti cle porosity abets stickin g by allowing for greater dissipation of kinetic 
energy (jPominik & Tielens 1997h . Moreover porosity lowers p s , which permits 



larger masses to coagulate when growth is limited to a fixed r s , as above. The 
degree of enhanced mass growth depends on the drag regime; fixing the stopping 
time gives a particle mass m oc p~ 2 for Epstein drag, and m oc 1/^/Ps for Stokes 
drag. These advantages are limited, however, since collisions at speeds of 0.1— 
lms" 1 restructure and co mpactify aggregates composed of //m-sized monomers 
(iDominik fc Tielens 19971 ). 



Simulations of grain growth by Dullemond fe Dominik f2005h and Ormel et al. (20071 ) 



include vertical settling and neglect radial drifts. iDullemond &: Dominik (2005) 



find that in laminar disks, porous aggregates achieve equivalent compact sizes 
of a few cm before settling to the midplane, while particles assumed to be al- 
ways compact gro w up to 1 cm (see their Figure 3, models S2, S5, and S6). 
Ormel et al. (20071 ) calculate that in turbulent disks, porous bodies attain equiv- 



alent compact sizes of several cm, for a = 10 2 at r = 1 and 5 AU. 

In both these studies, /um-sized grains deplete, and the disk becomes optically 
thin, within ~10 3 yr. This timescale is too short to be reconciled with astronom- 
ical observations of disk spectra. To maintain the population of small grains and 
the disk's optical depth over 1-10 Myr, models may need to account for collisional 
fragmentation of ag gregates, or condensation of dust from silicate v apor in hot, 
active disk regions (jDullemond fc Dominik 20051 . lOrmel et al. 20071 ). Turbulent 



concentration of dust ( §3.2.3p could also be an important correction to coagula- 
tion models, not only because densities are enhanced, but because collision speeds 
are likely r educed in particle clum ps, as measured in simulations of the streaming 
instability (j.Tohansen et al. 200 9b. g72]). 



5 CRITERIA FOR GRAVITATIONAL INSTABILITY 
5.1 Dynamical Collapse 

We analyze the stabili ty of a planar, s elf-gravitating, rotating sheet of dust, mod- 
eled fluid (GW, ISafronov 1969h . For the moment we ignore interactions 



with gas. Unperturbed, the dust has surface density E, angular speed O, and 
barotropic pressure -P(S) ~ £c 2 , where c is the velocity dispersion. For axisym- 
metric perturbati ons, the linearized equation s for continuity and momentum read 
(e.g., chapter 6 of Binney &: Tremaine 20081 ) 
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where perturbations are primed, $' is the perturbation gravitational potential, 
and k = (rdVl? jdr + 4f2 2 ) 1 / 2 is the epicyclic frequency of radial oscillations — in 
practice nearly equal to f2, its value for a Kepler potential. Fo r WKB waves, per- 



turbations oc 

e t{k r r+ut) w here£; r r > 1, and « -2vrG£7| fc r | dBinnev & Tremaine 200 



Then (|18p -(|20 p yield the well-known dispersion relation for axisymmetric waves: 

uj 2 = c 2 k 2 r -2ttG£|Av| + k 2 (21) 

which informs us that pressure stabilizes short wavelengths, rotation stabilizes 
long ones, and self-gravity de-stabilizes intermediate wavelengths when 

Q = < 1 (22) 



(jToomre 19641 . iGoldreich Lvnden-Bell 19651 ) . Modes whose wavelengths ex- 



ceed A cr j t = Ai{ 2 GYi/k 2 (neutrally stable for c = 0) are always stable. The fastest 
growing mode has Af gm = 2c 2 /(G£) = Q 2 A cr it/2. 

A standard expression for the vertical thickness of the dust layer is h p ~ c/il 
(caveat emptor: (|16p shows this is not valid when St <C 1). Then Toomre's 
criterion (|22p translates into a thickness criterion 

ttCY f r \ 3 / 2 

^ < K « ^ « 2 x 108i?Z -i (^j cm (23) 
or equivalently, a requirement on the midplane density 



-3 



'^"^"S" 18 "'^ gcm " (24> 

which has the same scaling as the Roche criterion (31]), but is smaller by a geo- 
metric factor of ~20 that arises because the gravity of an axisymmetric ring is 
greater than that of a tidally distorted ellipsoid. The conditions for dynamical 
collapse differ by orders of magnitude from those afforded by the gas disk: 

/ r \ 3 / 14 

h;/h g « 5 x W- 4 FZ rcl J (25) 

35 / r \~ 3 / 14 
,7P S - f (^) • (26) 

Turbulence would have to be exceptionally weak to allow particles to sediment to 
such a thin layer: from ((HI), a would need to be ^ 3 x 10 _7 r s (FZ re i) 2 (r/AU) 3/14 . 
Overcoming this obstacle — and showing that there are alternatives to dynamical 
collapse from a disk of uniform Solar abundance — are the goals of this review, 
with some of the main ideas summarized in Figures [2] and [4l 

Sekiva (1 983) solves for the linear stability of a midplane layer composed of 



perfectly coupled dust and gas. The layer is confined by the pressure of over- 
lying gas layers. As with Toomre's calculation above, axisymmetry is assumed; 
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however the midplane layer's vertical thickness is free to change. Dynamical 
instability occurs for midplane densities about 4 times greater than p* because 
gas pressure helps to stabilize the midplane layer. Since the layer is incompress- 
ible, in-plane motions give rise to out-of-plane "bulges." Though not modeled, 
slow sedimentation of dust to the centers of these bulges could produce km-scale 
planetesimals. 



Toomre's criterion (|22p is a good rule of thumb for deciding when self-gravity 
matters, but it misleads because it is derived for axisymmetric waves. In reality, 
GI proceeds non-axisymmetrically. Non-axisymmetric waves amplify when swung 
from leading to trailing by the background radial shear ( Goldreich &; Lvnden-Bell 19651 ). 



nvcs 



Waves grow even if Q > 1, but especially strongly as Q — > 1. iToomre (198ll ) g 
a delightful tutorial, explaining swing amplification as a near-resonance between 
radial shear, epicyclic motion, and self-gravity. Amplification is restricted to 
the interval, of duration ~l/fJ, when wave pitch angles slew from about -1 to 1 
rad; upon completion, A r ~ and amplitudes have grown by factors ^ 50 for 
Q 1.2. Modes most prone to growth have « A cr it- For these reasons, GW 
describe the disk's initial fragments as having a characteristic length Lf rag ~ £A cr it 
and mass Mf rag ~ ££ 2 A 2 rit = £ 2 M cr i t , where £ is an order unity parameter that 
contains our uncertainty about the spectrum of seed perturbations and how close 
Q is to unity. At r = 1 AU, Mf rag equates to a compact rocky planetesimal having 
size Sf rag 8(FZ re i/0.33)£ 2 / 3 km; this estimate figures prominently in the lore of 
"kilometer-sized planetesimals." 

Numerical N-body simulations of inelastically colliding par ticles verify that 
within 1-2 orbital periods, particles aggregate on the scale A cr it ( Michikoshi et al. 20071 ). 



The li near mass scal e M CI [ t is not as prominent in either particle or gas simula- 



tions (jGammie 20011 ) , partly because fragments accrete rapidly — suggesting that 



-Merit be used cautiously, and probably as a lower limit. 



5.2 Drag- Assisted Gravitational Instability 

The stability properties discussed above pertain to dust treated as a single fric- 
tionless fluid. But gas-dust interactions can change this picture qualitatively. 
The s implest modificat i on is to introduce d r ag terms to the momentum equa- 
tions dWard 197fil . l2000l . iGoradini et al. 198ll . lYoudin 2005h : 



dv' -c 2 dXl' d& v' 

^ - = —— - — -Hl (27) 
8t * £ 8r dr t B { 1 

The meaning of v 1 bears clarification. With drag, particles drift radially ( ^3.1 . 1 j) . 
Therefore perturbed quantities refer to a background with steady v r ^ 0, and 
our WKB analysis is restricted to t < t r . Backreaction on gas, whose properties 
are assumed fixed, is neglected. 

When t s <1, the drag-modified equations yield a new dispersion relation: 

uj = i(c 2 k^. -2irGT l \k r \)t s (29) 

which implies that modes for which k r < 2irGT,/c 2 are unstable — even when 
Q > l! The fastest growing mode has Af gm = 2c 2 /(GS) = Q 2 \ C rit/2 and growth 
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rate |wf gm | = Qt s /Q 2 . Thus drag destabilizes long wavelength modes, on longer 
timescales, compared to the dissipationless case. These qualitative differences — 
in particular the removal of Toomre's Q-criterion — apply for a variety of drag 
regimes, and for cases that include turbulent diffusion of particles, accounted for 
by adding (Z?p/£)<9 2 £'/<9i 2 to the right-hand side of (PT8j) (Youdin, in preparation; 
Karim Shariff & Jeffrey Cuzzi, personal communication). 



Goodman &: Pindor (20001 . hereafter GP00) explain the instability simply. In 
the unperturbed state, dust has constant surface density and drifts radially in- 
ward at constant velocity. An axisymmetric overdense ring of width Ar exerts 
a gravitational pull radially inward at its outer edge. Provided h p <C Ar <C h g , 
dust rotates faster there, while the gas velocity is unaltered. The increased drag 
causes dust to flow at a greater rate into the annulus. Likewise there is a flow 
into the annulus at the inner edge, where dust experiences less drag. Thus when 
combined with gas drag, self-gravity, no matter how weak, can draw particles 
into rings. 

Th e stability of non -axisymmetric modes with dissipation is a largely open 
field. iNoh et al. Q99lh make some exploratory integrations of the linearized per- 



turbation equations. In N-body simulations that include —v'/t s drag and a fixed 
gas velocity field, particles cluster m ore readily because gas drag damps their 
velocity dispersion ( Tanga et al. 2004] ) . 



6 THIN DUST LAYERS: VERTICAL SHEARING 
INSTABILITIES 

Particles may be prevented from settling to the midplane not only by turbulence 
intrinsic to the gas (e.g., su stained by the MRI) , b ut also by turbulence triggered 
by the particles themselves ( Weidenschilling 1980l ). We consider here instabilities 



caused by small r s <C 1 solids that drive a vertical shear. In $7] we show that 
larger particles, in addition to inducing a vertical shear, can drive turbulence of 
a different character via relative streaming motions. 

As dust settles, a particle-rich sublayer develops which orbits at nearly the 
full Keplerian velocity. The dust-poor gas above and below shears by at speeds 
approaching the pressure-supported value tjvk (eq. E|). Such a stratified shear 
flow is subject to a Kelvin-Helmholtz-type instability (KHI). For a non-rotating 
flow — which we later show turns out to be a fair analogue for actual vertically 
shearing flows in disks — the KHI can develop when the Richardson number 



somewhere (jDrazin Sz Reid 20041). where q is the vertical gravitational acceler 



ation and p = p p + p g (see iGaraud Lin 20041 . for the applicability of Ri to 



a dust-gas mixture). The square root of Ri is the Brunt- Vaisala frequency for 
buoyant vertical oscillations, divided by the shear rate dv^/dz. This is a sen- 
sible criterion: when the buoyancy frequency is small, vertically displaced fluid 
elements fail to return to their equilibrium positions before being whisked away 
by the background shear. 

When we neglect self- gravity and approximate g ~ —Q^h p , and further take 
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d{\xip)/dz ~ — l//i p and \dv ( j ) /dz\ ~ r/v^/hp, criterion ([30]) implies that the KHI 
prevents the sublayer from having a thickness less than 

/i p , Ri ~ Ri 1 / V ~ lO" 2 (^j) ^ >»g • (31) 

Thus the midplane densit y falls short of p* b y a factor of ~20/(i ? Z re i), roughly 
consistent with findings by Cuzzi et al. (19931 ) , who use mixing length models for 



particle-laden turbulence scaled to laboratory experiments. If we insist that the 
disk be of bulk Solar composition, then achieving higher densities requires either 
that particles be larger (r s ^ 1) so as to decouple from the turbulence, or that 
they clump locally f 33X31 g7|). 

Nominally we can estimate, using (|16p . that turbulence intrinsic to gas dom- 
inates the KHI when a £ T s (/i PiRi //j g ) 2 ~ 3 x Kr 4 T s (r/AU) 4 / 7 . Whether MRI- 
dead zones ($2]) are sufficiently passive for the KHI to manifest remains an out- 
standing question. Our a scalings may not apply to the weak motions possi- 
bly excited by M RI-activ e surface layers ( Gammie 19961 . Fleming &; Stone 2003 . 



Oishi et al. 20071 : see also iBai Goodman 20091 for detailed estimates of the ac- 



tive layer thickness). 

In one of the earliest works to t reat the KHI in disks re alistically by account- 
ing for orbital differential rotation, llshitsu &: Sekiya f2003h numerically integrate 



the linearized perturbation equations and find that maximum growth factors in- 
crease strongly with midplane dust-to-gas ratios — so strongly that p p / p g is likely 
to stall at values just above unity in disks having b ulk solar meta l licity. Their con - 
clusion is borne out by shearing box simulations (jChiang 200 «1. iBarranco 20091 ). 



performed in the limit that dust is perfectly coupled to gas (r s — > 0). These 
simulations demonstrate that despite the inclusion of orbital dynamics, Richard- 
son numbers characterizing instability remain between 0.1 and 1 (Figure [5]). In 
retrospect this is not surprising, since rotational and radial shearing frequencies 
~J1k are only of the same order as both the Brunt- Vaisala frequency, and the 
vertical shearing frequency h p ^i/{r]V%) ~ Ri 1 / 2 !^ 1 , when Ri ~ 1. Nevertheless 
radial shear is critical for limiting the growth of p erturbations that otherwi se run 
away when only the Coriolis force is included (cf. Gomez &: Ostriker 20051 ). 



As pointed out by J. Goodman (personal communication; see Chiang 20081 ). 
the KHI may not be the only instability at play. Dusty sublayers may also be 
baroclinically unstable: their isodensit y surfaces, dominated by dust, do not align 



with their isobars, dictated by gas (Knobloch &: Spruit 19851 s ) . Axisymmetric 



baroclinic stability is assured by Ishitsu &: Sekiva ('20031') and by analogy with the 



Solberg-Hoiland criteria for rotating stars ( Kippenhahn h Weigert 1990l ). Specif- 
ically, displacements at constant entropy are almost perfectly radial when vertical 
stratification of the dust is strong. Thus in a Keplerian disk, outward adiabatic 
displacements move toward regions of higher specific angular momentum, and 
are Rayleigh stable. 

Aside from postulating large, aero dynamically decoupled particles, another way 
to p revent the she aring instability from forestalling GI is to substantially raise 



J re\- 



Sekiva (19981 ) discovers that for flows with constant Ri(z) = 1/4, Z xc \ rs 6- 



30 can lead to p(z = 0) rs p* in a midplane cusp. In fact, by further increasing 
Z Te \ slightly, and accounting properly for vertical self-gravity, the cusp becomes 
singular: p(0) — > oo. YS perform similar calculations, finding p ~ p* could be 
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achieved for Z rc \ as low as 2, depending on disk parameters. They show that 
the cusp develops when p p ;> p g throughout the layer, which from (|3ip requires 

Ep/Eg ^ \JW\tivk/ c g , or Z TC \ ^ l(r/AU) 2 / 7 . Physically, as solids are added, the 
buoyancy frequency increases, while the shear saturates at rjVK/hp. Numerical 
simulations, which as yet do no t include vert ical self-gravity, suggest Z TC \ ;> 5 to 
attain densities approaching p* ( Chiang 20081 ). Mechanisms to achieve supersolar 



metallicities are shown in Figure [2] and described in §3.1.21 

Whether or not they are gravitationally unstable, thin dust layers for which 
p p ;> p g are also prone to the streaming instability ( §7.2p . a mechanism that can 
further concentrate solids relative to gas when t s / 0. 



7 DRAG INSTABILITIES 

The vertical shearing instability described in £j6] arises in dense particle layers, 
even when dust is perfectly coupled to gas (r s — ► 0). When relative motion 
between gas and dust is allowed (r s ^ 0), new secular (i.e. requiring dissipation) 
instabilities arise. These drag instabilities can strongly clump particles relative 
to gas. As with vertical shear, drag instabilities require that p p /p g be large 
enough that particles backreact significantly on gas. Self-gravity is not required, 
in contradistinction to the drag-assisted GI of §5.21 



We describe two different models of drag instabilities. The secular dust layer in- 
stability (SDLI) of GP00 considers the drag arising fr om turbulent stresses on the 
stratified dust layer. The streaming instability (SI) of Youdin &, Goodman (20051 . 



hereafter YG05) accounts for the drag forces acting mutually between gas and 
particles in a laminar unstratified disk. Both analyt ic models — as well as numer- 
ical simulations that add extra physics to the SI ( Johansen et al. 20071 ) — show 
that particles clump. 



7.1 Secular Instability of the Dust Layer 

The SDLI uses a single-fluid, height-integrated, axisymmetric model of the dust- 
rich sublayer. Dust-poor gas at the top and bottom of the sublayer shears by 
at relative speed tjvk, and is envisioned to be turbulent ($6]). These turbulent 
boundary layers exert drag on the dust layer. Turbulent surface stresses are 
assumed to be communicated throughout the entire layer; the detailed vertical 
response is, by construction, igno red. To calculate the turbulent drag, GP00 



adopt prescriptions from GW and lCuzzi et al. (19931 ) for "plate drag," drawing 



an analogy with Ekman flow past a rigid plate. The SDLI generates overdense 
rings on orbital timescales, with radial widths comparable to the sublayer vertical 
thickness. 

The use of the plate drag form ula is not well founded, e xcept perhaps if particles 



in the sublayer all had r s ~ 1 (jYoudin Chiang 2004 ). But GP00 argue that 



their results are more general, transcending the specific drag prescription used 
to obtain them. The crucial assumption underlying the SDLI is that drag is 
collective, i.e. it depends on the surface density of the sublayer. GP00 construct 
a toy model that demonstrates how any collective drag force leads to instability. 
The insights apply to drag instabilities generally, including the SI. 
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In the toy model, the sublayer is described by a ID distribution of mass density 
£ and velocity v: 

(32) 

g-v d {Y,)v (33) 

where fd(£) > accounts for drag, and g < serves as a proxy for gravity, 
pressure, and Coriolis forces. The equilibrium state has uniform mass density So 
and inward drift velocity vq = g/u d (T, ) = gjv§. We might imagine v& increases 
with £ — insofar as the plate drag torque is independent of £ while the angular 
momentum of the sublayer is proportional to £, so that the drift speed decreases 
with increasing £ — but the details are unimportant: instability arises as long 
as v& and g depend on £ differently. With no loss of generality, all of the £- 
dependence is relegated to v&. 

We take linear perturbations, denoted by primes, to have a Fourier dependence 
exp(rt — ikx). The drag coefficient is expanded 

z/ d (£) = ^ d (£ ) + dv d /dZ\ £' . (34) 

A non-zero drift vq is necessary for perturbations to u d to enter linearly. The 
dispersion relation reads 

r = |f-l + 2i~k ± ^1-Ai5 u kj , (35) 

where k = kv^/vQ and 5 U = din fd/dln £| . For any 5 U ^ the positive root gives 
growth, since K ^\/l + > 1 for all real b. Thus any collective drag produces 
instability. 

To understand the growth mechanism, we examine the eigenfunctions £' and 
v', first Taylor expanding the growing mode about 5 U , 

T = T + Ti + T 2 = ikv - i5 u kv + v k 2 5l , (36) 

and similarly expanding v' = v'q+v^+v^- By substituting these into the linearized 
eqs. ([32]) ~ (j33|) and equating terms of the same order, we see that v' = 0, v[ = 
— 5 I/ (£'/£o)fo, and v' 2 = —(Ti/vd)v[ = — «M^(£'/£o)wo- For 5 U = 0, oscillations 
are stable with wave speed ^s(To)/k = vo. This is a neutral mode where £' simply 
advects with the background flow. To first order in 5 U , v[ is anti-phased with 
£'; if 5 U > 0, the drift speed decreases (increases) at density maxima (minima) 
due to the perturbed drag force; the wave does not amplify but its speed shifts 
by Q(Ti)/k = —5 u vq. To second order in 5 U , the acceleration of v[ is subject to 
drag, which induces a secondary flow v' 2 . Since v' 2 is — tt/2 out of phase with £', 
the mode grows, at rate T 2 - The phase shift enters because accelerations change 
sign at density extrema. Figure [6] depicts the eigenfunctions. 

The toy model shows that essential components of drag instabilities include the 
background drift vo, which allows collective effects to enter linearly, and time- 
dependent oscillations that drag forces overstabilize. 



d£ d .„ , 



dt 



dv dv 
dt dx 
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7.2 Streaming Instability 

Since the streaming instability involves a minimal amount of physics — Keplerian 
orbital motion, gas pressure, and a drag acceleration that is linear in relative 
velocity — it supports the robustness of drag instabilities. YG05 uncover the SI 
by modeling gas and dust as two interacting fluids that obey 

= -^-^T 1 (37) 



Dt * p g t 



s 



Dj : -PpV-^p (39) 
V-v g = (40) 
where Dj / Dt = d/dt + Vj-'V. Vert ical gravity is neglected, so r is the cylindrical 



radius. lYoudin Johansen (20071 . hereafter YJ07) confirm that gas compress- 
ibility is justifiably neglected in (Ti0j) . and that the fluid approximation for dust 
holds! 

The steady-state solutions to (f3"7"|) - (l4T)l) give Keplerian moti on, plus radial and 



azimuthal drifts between dust and gas (INakagawa et al. 1986 . §3.1. lj) . The rel- 



ative streaming motion is linearly overstable (YG05), behaving similarly to the 
GP00 toy model. The main complication is that while unstable modes can be 
axisymmetric, they necessarily involve motions in all three directions. 

Particle clumping occurs through the backreaction term in (|38|): —p p {v p — 
v g )/(p g t s ). Since p g and t s are constant — by assumption in YG05, and to excel- 
lent approximation for subsonic flows — only perturbations to p p couple linearly 
to background streaming motions v p — v g . Thus particle clumping is required 
to extract energy from streaming motions, as in the GP00 toy model. The SI is 
ultimately powered by the background gas pressure gradient, which does work on 
radially flowing gas: —v gr dP/dr > 0. Since streaming motions transport angular 
momentum inward, no energy is extracted from the Keplerian shear flow (YJ07). 

The SI is controlled by two parameters: r s and (p P /pg)o, the background dust- 
to-gas ratio. While the SI always exists, growth is slow in the test-particle limit 
(p p /pg) -C 1, and in the perfect coupling limit r s < 1. As r s — ► 1, background 
drift speeds peak, and modes can grow on orbital times. 

The nonlinear consequences of the SI are best studied by computer simula- 
tions, which can selectively include other physics. Treating particles as a con- 
tinuous fluid simplifies analytic calculations, but causes difficulties in numerical 
simulations. The use of artificial viscosity to avoid density discontinuities in 
the pressureless fluid can underestimate particle concentrations and otherwise 
compromise results. Alternatively, "hybrid" simulations model g fluid on 

an Eulerian grid and solids with Lagrangian superparticles, each representing a 
swarm of actual particles. YJ07 simulate the linear growth of SI, and confirm 



3 We are not aware of a formal criterion for deciding when dust can be modeled as a fluid. 
The usual criteria for gas molecules, that they be collisional, do not apply to dust particles 
entrained in gas. YG05 suggest that coupling to the gas may suffice: if ujt s <C 1, then the fluid 
approximation is expected to be valid on ti mescales a? -1 . A more rigorous criterion would be 
useful since many nonlinear simulations (e.g. Iliogan fc Cuzzi 20071 ) use the two-fluid approach. 
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that the hybrid and two-fluid approaches give convergent results. They also pro- 
vide the eigenfunctio ns of growing Four i er modes, useful fo r testing codes with 
two-way drag forces ( Balsara et al. 20091 . Ilshitsu et al. 20091 ) . 



Johansen Sz Youdin (20071 ) perform nonlinear 3D simulations of the idealized 
SI using superparticles. Marginally coupled, r s = 1 solids concentrate by factors 
of several hundred for a range of 0.2 < {p p /p g ) < 3.0. Better coupled r s = 0.1 
particles give particle overdensities of several tens, but only for {p p /p g ) > 1. 
When {pp/pg) = 0.2, particle clumping is of order unity. Sample simulations are 
displayed in Figure [HJ 

Realistically assessing how effectively the SI can concentrate par t icles r equires 



that vertical stratification be included (e.g. iJohansen et al. 2006al . 120071 . l2009bl . 
hereafter JYM09), since high dust-to-gas ratios also trigger the KHI ($6]). JYM09 
perform 3D hybrid simulations of r s = 0.1-0.4 solids to study the combined 
effects of KHI and SI. They find a metallicity threshold for strong clumping: 
when Z re \ increases from 2/3 to 4/3, the SI, unaided by self-gravity, triggers 
overdensities having p p > 10 3 p g . The Z re \ threshold corresponds to (p p /p g ) ;> 1 
at the midplane, similar to the analytic criterion for the saturation of the KHI 
(YS). It remains to be seen if this threshold will shift, or be manifested differently, 
for smaller particles. Hybrid simulations with SI that include self-gravity are 
discussed in 351 



8 GRAVITATIONAL COLLAPSE INTO THE NONLINEAR 
REGIME 

The fragmentation of the disk into marginally self-gravitating clumps is only the 
beginning of the story of planetesimal formation. About six orders of magnitude 
need to be traversed from the Roche density to solid densit y. Even if clump s 



survive buffeting by turbulence and ram pressure stripping (jCuzzi et al. 20081 ). 
their further contraction is resisted by internal random motions and net angular 
momentum. Gas drag slows the sedimentation of particles toward clump centers, 
but it also enables such settling by braking the clump's rotation and damping 
random motions. For t s < 1, collapse proceeds on timescales longer than the 
free-fall time; how much longer is unclear. At large enough densities, inelastic 
collisions between particles further reduce random motions and accelerate gravi- 
tational collapse. 

The popular notion that GI produces planetesimals having sizes on the order of 
a kilometer probably originates from GW, who recognize that a fragment having 
the full mass Mf rag ($5]) has too much angular momentum to collapse unim- 
peded to solid density. They propose instead that first-generation planetesimals 
form from overdense regions having horizontal sizes X p r 9 / 4 /(/5s^ 4 M* -C A cr j t ; 
such clumps, if they conserve angular momentum, collapse into solid objects 
having sizes S ~ QA FZ rP } km, spinning just below break-up (see also §15.2 of 
Goldreich et al. 2004h . But this estimate is unjustified because it assumes that 



scales <C A cr j t can collapse; this would require Q <C 1, and indeed GW assume 
unrealistically that c = 0. 

While reliable analytic estimates of the sizes of first-generation planetesimals 
are lacking, numerical simulations have made great strides, providing a proof of 
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principle that gravitational collapse can occur, even in turbulent flows. iJohansen et al. (20071 ) 



synthesize practically all of the physics described in this review by executing 3D, 
self-gravitating, vertically stratified, shearing box simulations of superparticles in 
MRI-turbulent gas. Particles all have r s = 0.25-1. Three ways of concentrating 
particles manifest: trapping within transient maxima in gas pressure ( §3.1. 2ft ; 
the streaming instability ( §7. 2D : and self-gravity. Within only a few orbital peri- 
ods, the largest bound cluster of superparticles attains a mass ~3x that of the 
450-km-radius asteroid Ceres — though collapse to solid densities is not explic- 
itly followed because of finite grid resolution. Surprisingly, magneto-rotational 
turbulence hastens gravitational collapse, through pressure trapping. The main 
uncertainties of these simulations are the assumption that sticking produces the 
largish seed particles, having sizes of 15-60 cm at 5 AU; and the specific real- 
ization of MRI turbulence, which depends on the unknown net magnetic flux 
and on microscopic dissipation parameters that are too small to si mulate exactly 
(jLesur & Longaretti 20071 . iFromang et al. 20091 . iDavis et al. 20091 ) . 



MRI turbulence, though it appears to aid gravitational collapse, is not a nec- 
essary ingredient. Figure [8] displays how bound clumps also form in an unmag- 
netized simulation (JYM09). 



9 SUMMARY AND OUTLOOK 

For a collection of grains within a protoplanetary disk to become gravitationally 
bound, their density p p must exceed the Roche density, which for typical param- 
eters is over 1000 times the ambient gas density p g . Equivalently, the dust-to- gas 
ratio must be enhanced over the protosolar value by a factor of 60,000. Daunting 
as it may seem, the requirement may be met by several mechanisms, working in 
concert, that concentrate particles relative to gas. 

It certainly helps to have a disk — or regions of the disk — with an abundance 
of solids relative to gas greater than that given by solar composition. The dust- 
to-gas colum n density ratio S p /Sp- may need to be supersolar by factors of 2-10 



(£j6j see also IJohansen et al. 20071 ) for planetesimals to form. The sharp ris e in 



exoplanet detection rates with host star metallicity ( Fischer Valenti 20051 ) re- 



inforces the importance of initial dust abundance. Disks can also become enriched 
with time, as pile-ups due to radial drifts and photoevaporation of gas increase 
Sp/Sg (Figure [2|). Observations of gas giants are consistent with their origin in 
metal-enriched disks. B ulk metallicities range from ~4x solar (Jupiter; JQ) to 
~50x solar (HD 141569: ISato et al. 20051 ). 



Two popular misconceptions deserve correction. First, self-gravity need not 
be comparable to stellar tides before its effects are felt. Gas drag on particles, 
abetted by arbitrarily weak self-gravity, collects particles into overdense axisym- 
metric rings ($5]). Second, gas turbulence does not necessarily obstruct GI. There 
is growing evidence that magneto-rotational turbulence creates long-lived pres- 
sure maxima that can trap particles ( §3.1.21 §3.2j) . 

The above concentration and enrichment mechanisms — combined with vertical 
settling (Figure [2]) — can eventually lead to p p /p g ~ 1. This condition presents 
a much lower hurdle than the Roche criterion, or even Toomre's Q-criterion. 
Having as much inertia in particles as in gas paves the way for collective drag 
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effects — notably the streaming instability — that can clump particles still further, 
all the way to Roche density. Figure H] illustrates the size ladder that dust must 
climb to form planets, identifying key processes. 

Below we chart possible courses for future work: 



1 . Chondr ules — arguably the building blocks of first-generation planetesimals — 
were flash-heated in an extremely dust-rich environment, with p p possibly 
exceeding p* (SID- That none of the concentration mechanisms discussed 
in this review (see §3.2.30 explains how they were heated indicates our 
story may be incomplete. In shock models for chondrules, the isothermal 
post-shock gas density exceeds the pre-shock density by the square of the 
shock Mach number, a factor of about 30 — just enough for p g ~ p*. This 
opens the exciting possibility of a violent origin for both chondrules and 
planetesimals, though leaves unidentified the source of such strong shocks. 

2. To what extent does dust settle to the disk midplane? The answer depends 
on stirring by turbulence, whether driven by the solids or externally im- 
posed. Common idealizations — notably that the dust density profile in oth- 
erwise laminar disks has constant Richardson number Ri ~ 1/4 — hold only 
approximately at best, as shown by s imulations pe r formed in the li mit that 
dust and gas are perfectly coupled ( Chiang 2008 . Barranco 20091 ). Long- 
term evolution of the dust distribution requires relaxing the perfect cou- 
pling assumption, so that dust can sediment out of gas. Novel numerical 
approaches such as implicit time evolution and multiple time stepping will 
need development. The wide separation between stopping, orbital, and sed- 
imentation timescales could reveal new phenomena. There may not even be 
a steady state to which dust settles. In laminar disks, the dust distribution 
may cycle through long periods of settling interrupted by brief bursts of 
nonuniform turbulence wherever Ri drops below unity. 

3. Simulations that include the streaming instability (e.g., Johansen et al. 20071 ) 
find that particles moderately coupled to gas — having stopping times r s > 
0.1 — clump strongly enough to undergo gravitational collapse. Does the 
same conclusion apply to yet smaller solids, increasing the overlap between 
process 1 — growth by sticking — and process 3 in Figure UP Pessimistically, 
smaller solids might not settle to a layer with p p /pg <; 1 where streaming in- 
stabilities are most effective. Also clump lifetimes decrease with decreasing 
stopping time (Figure [71 §3.2.3p . Optimistically, particle concentration is 
observed in current simulations to increase with finer grid resolution, with 
no convergence as yet. An adaptive mesh applied to a shearing box would 
be a powerful tool to track the densest particle clumps. 



Together points 2 and 3 emphasize that elucidating the respective roles of 
sticking and self-gravity requires bridging our understanding of particle dy- 
namics from the aerodynamically well-coupled to the marginally decoupled. 

4. The post-Roche evolution of particle clusters is largely unexplored. A clus- 
ter could fragment into multiple ones, much as interstellar gas clouds un- 
dergo hierarchical fragmentation as the Jeans mass decreases with higher 
density. The number of fragments could be large if the initial overdensity 
is an azimuthally extended ringlet. A record of the initial mass function for 
planetesimals might even be preserved today, in the mass functions of as- 
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teroids and Kuiper belt objects dMorbidelh et al. 20091 ). Extremely widely 
separated binaries ( Petit et al. 20081 ) may be gravitationally collapsed clus- 
ters that fissioned. 
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Figure 1: Hubble Space Telescope image of t he protoplanetary di sk HH 30, seen 
in scattered starlight at optical wavelengths ( Burrows et al. 19961 ). Overlaid are 



contours from 1.3- mm continuum radiatio n measured by the Plateau de Bure In- 
teferometer (PdBI. lGuilloteau et al. 20081 ). The optical image reveals the opaque 
disk midplane, a dark lane dividing two "bowls" of light scattered by sub-/xm 
grains in flared disk surfaces. Perpendicular to the disk plane is the jet. Mm- 
wave emission arises predominantly from the midplane, plausibly from cm-sized 
grains that have settled vertically. The mm- wave disk has an outer radius of ~130 
AU, smaller than the optical radius of ~300 AU, suggesting either that grains 
grow faster at smaller radius or that cm-sized grains have drifted radially inward 
( §3.1.1 j) . Axes are measured in arcseconds; 1" = 140 AU for an assumed distance 
of 140 parsecs. The inset shaded spot shows the size of the PdBI beam, a mea- 
sure of the angular resolution of the mm-wave map. Reproduced by permission 
of Astronomy & Astrophysics. 
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Figure 2: Four mechanisms for metal enrichment — none of whic h involve self- 
gravity. Dust can settle vertically into a dense sublayer (e.g., Chiang 2008, 
Barranco 2009. fl6l): pile up as it drifts radially (jYoudin h Shu 20021 . 
Youdin Chiang 2004 . ^3.1.21); remain behind as stellar ultraviolet radiation 



photoevaporates gas (e.g., Throop &: Bally 20051 ); and be concentrated on small 
scales by passively responding to turbulent fluctua tions (£13. 1.21 ^3.2.31) and 
by actively driving drag i nstabilities with gas (e.g ., 
Youdin Goodman 20051 . Ijohansen & Youdin 2007I . 



Goodman fc Pindor 20001 . 
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Figure 3: Test particle s concentrate in turbule nt pressure maxima, from 3D shear- 
ing box simulations of ljohansen et al~ The surface densities of gas (left) 
and r s = 1.0 solids (right) are azimuthally-averaged and plotted versus radius x 
(in units of h g ) and time t (in orbits). Magneto-rotational turbulence drives gas 
fluctuations, to which particles passively respond via drag forces. The order-unity 
particle density fluctuations overlap the ~l%-scale gas density fluctuations — 
equivalent to the pressure fluctuations for the assumed isothermal gas. Careful 
inspection shows that solids collect slightly downstream (smaller x) from local 
pressure maxima since this is where local gradients dP/dx > cancel the global 
dP/dr < (see §3.1.2h . Reproduced by permission of the American Astronomical 
Society. 
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Figure 4: Scaling the size ladder from dust to planets. Physical processes that 
either grow (labeled with white circles) or concentrate solids (yellow) are illus- 
trated. Numbers advance from the earliest to latest stages of planet formation. 
The relevance of a given mechanism tends to be restricted to a certain range of 
particle sizes, indicated crudely by the arrows on the ladder. The ranges shown 
are subject to debate and actively researched. The least well explored is shown by 
a dashed arrow: drag-assisted GI ( §5.2|) . Dimensionless stopping times r s = Or^s 
are shown for r = 1 AU; some processes are aerodynamic and depend more on r s 
than on particle size. 



PLANETESIMAL FORMATION—SUBMITTED DRAFT 



35 



N 



CD 




-0.01 



-0.04 



-0.02 



y(h g ) 



Figure 5: Dusty sublayer undergoing a Kelvin-Helmholtz-like instability, taken 
from a 3D shearing box simulation by Lee et al. (in prep.), using the anelastic 
code of Barranco The Richardson number is set initially to a constant 

Ri = 0.1 ( Sekiva 19981 ). with a midplane dust-to-gas ratio of 10. Larger initial 
values for Ri ^ 1 produce no instability for at least 10 orbits. The vertical 
shearing parameter rjVK/cg = 0.025. The size of the shearing box is (L x , L y , L z ) = 
(12.8,6.35,8) 

^max) where z max is the maximum height of the dust layer (e.g., 
Chiang 20081 ) . and the number of grid points is (N x , N y , N z ) = (128, 32, 128). 
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Figure 6: Me chanism behind particle c lumping in drag instabilities, using the 
toy model of Goodman &; Pindor (200ol ). The equilibrium inward drift speed vq 
(grey line) is constant with local radial coordinate x. For linear modes, a surface 
density perturbation X' (solid curve) of wavelength 2ir/k has a related velocity 
perturbation v'. To first order in 5 V = d\n i^/dln S > (eqs. [33Tl34p . overdense 
regions slow their drift (v[, dotted curve). Since resultant mass fluxes (dotted 
arrows) point toward density zeroes, there is no amplification of £' at this order; 
the mode merely oscillates. The second-order response v' 2 (dashed curve) arises 
from drag forces acting on the first-order flow, and is phase-shifted relative to 
There are now mass fluxes (dashed arrows) that amplify the original density 
perturbation; the mode is overstable. 
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Figure 7: Saturated states of the streaming instab ility, from unstratified sim - 
ulations of mutually coupled particles and gas by Johansen Youdin (20071 ) . 
Shearing box axes are oriented with radius running to the right, azimuth to 
the left, and vertical straight up. For the simulation on the left, t s = 1 and 
(Pp)/(p g ) = 0.2; for the right, r s = 0.1 and (pp}/(p g ) = 1- Colors show p p at 
the sides of the box, with black corresponding to zero density and bright yellow 
corresponding to densities up to ~30/9 g (left) or ~10/? g (right). Both snapshots 
are taken tens of orbits after the SI has saturated, when turbulence is fully de- 
veloped. For r s = 1, clumps are long-lived, so that Keplerian shear smears them 
into rings. Clumps of more tightly coupled particles are shorter-lived. Higher 
resolution simulations, which are feasible in 2D, produce still stronger clumping. 
Reproduced by permission of the American Astronomical Society. 
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Figure 8: Particle column density S p , showing the formation of seven gravita- 
tionally bound clumps in a 3D, vertically stratified, she aring box simulation o f 
unmagnetized gas and superparticles with r s = 0.1-0.4 ( Johansen et al. 2009bl ). 
The box-averaged (S p )/(S g ) = 0.02. The x(y) axis is parallel to the radial (az- 
imuthal) direction, and measured in units of the gas scale height. The simulation 
first evolved for 40 orbits without self-gravity, during which particles settled to 
the midplane and triggered vertical shearing and streaming instabilities; strong 
clumping resulted. This snapshot is taken 5 orbits after self- gravity was turned 
on. The bound fragments contain ~20% of the total mass in solids; each has a 
mass comparable to a compact planetesimal having a size 100-200 km. 



